/* Routine to pick a massive star from the Chabrier (2003) IMF (table 1), which is 
   psi(log(m))=.158*exp(-(log10(m) - log 0.079)^2/(2*.69^2)) (m <= 1)
   psi(log(m))=m^-1.3 (1 < m < 120)

We limit to the interval m>0.08


*********************************/ 
 

//Weight for the power low part or the exponential one
#define EXPPART  0.883623     

double chabrier_imf(gsl_rng *rand) {

  
  double x, mmin, mmax, gamma, outmass, mzero, sigma;

  x = gsl_rng_uniform (rand);
  
  if (x <= EXPPART) {
    //If in the lof part
    mmin=log10(0.08);
    mmax=log10(1.);
    mzero=log10(0.079);
    sigma=0.69;  //this has to be the stdev, not sqrt(2)*stdev    
    outmass=rangauss(mzero,sigma,mmin,mmax,rand);
    //from logm to m 
    outmass=pow(10.,outmass);
  } else {
    //If in the Salpeter part
    mmin = 1.;
    mmax = 120.;
    gamma = -2.3;
    outmass=ranpowerlaw(mmin, mmax, gamma, rand);
  }

  return outmass;
}
